"""
    AnisotropicFiltering

Advanced anisotropic filtering algorithms for flow-dependent background error 
covariance modeling in GSI. This module implements sophisticated correlation 
modeling that adapts to atmospheric flow characteristics, providing improved 
background error structure in regions of strong wind gradients and complex flows.

Anisotropic filtering addresses limitations of isotropic approaches by:
- Aligning correlation ellipses with atmospheric flow patterns
- Adapting correlation length scales based on local flow characteristics  
- Incorporating stream function analysis for flow direction determination
- Supporting multi-directional filtering with variable anisotropy ratios
- Integrating with existing recursive filtering infrastructure

Mathematical Framework:
The anisotropic correlation function is modeled as:
```
C(r₁,r₂) = σ²(r₁,r₂) · ρ(||r₁-r₂||/L(r₁,r₂,ψ)) · exp(-||r₁-r₂||²/2L²(ψ))
```

Where:
- L(r₁,r₂,ψ) is the flow-dependent correlation length scale
- ψ is the stream function characterizing the flow
- The correlation ellipse orientation follows local flow direction
- Anisotropy ratio adapts to flow deformation characteristics

Key Applications:
- Flow-dependent background error covariance construction
- Hurricane and tropical cyclone error modeling
- Frontal zone correlation enhancement
- Complex terrain flow adaptation
- Multi-scale atmospheric feature representation
"""
module AnisotropicFiltering

using LinearAlgebra
using SparseArrays
using StaticArrays

export AnisotropicFilter, FlowFieldAnalysis, CorrelationEllipse
export StreamFunctionAnalysis, FlowDirectionField, AnisotropyConfiguration
export apply_anisotropic_filter, compute_flow_direction, compute_anisotropy_tensor
export create_anisotropic_operator, validate_anisotropic_response

"""
    AnisotropyConfiguration{T<:AbstractFloat}

Configuration parameters for anisotropic filtering operations.

# Fields
- `base_correlation_length::T`: Base isotropic correlation length scale
- `anisotropy_ratio_max::T`: Maximum anisotropy ratio (L_max/L_min)
- `flow_adaptation_factor::T`: Flow adaptation strength parameter
- `ellipse_rotation_smoothing::T`: Smoothing factor for ellipse rotation
- `stream_function_filtering::T`: Stream function smoothing parameter
- `boundary_treatment::Symbol`: Boundary condition treatment (:periodic, :zero_gradient)
"""
struct AnisotropyConfiguration{T<:AbstractFloat}
    base_correlation_length::T
    anisotropy_ratio_max::T
    flow_adaptation_factor::T
    ellipse_rotation_smoothing::T
    stream_function_filtering::T
    boundary_treatment::Symbol
    
    function AnisotropyConfiguration{T}(;
        base_correlation_length::T = T(250_000),    # 250 km base scale
        anisotropy_ratio_max::T = T(4.0),           # 4:1 max anisotropy
        flow_adaptation_factor::T = T(0.75),        # 75% flow adaptation
        ellipse_rotation_smoothing::T = T(0.1),     # 10% rotation smoothing
        stream_function_filtering::T = T(0.05),     # 5% stream function smoothing
        boundary_treatment::Symbol = :zero_gradient
    ) where T
        new{T}(base_correlation_length, anisotropy_ratio_max, flow_adaptation_factor,
               ellipse_rotation_smoothing, stream_function_filtering, boundary_treatment)
    end
end

AnisotropyConfiguration(args...; kwargs...) = AnisotropyConfiguration{Float64}(args...; kwargs...)

"""
    CorrelationEllipse{T<:AbstractFloat}

Represents an anisotropic correlation ellipse at a grid point.

# Fields
- `center::Tuple{T,T}`: Ellipse center coordinates (x, y)
- `major_axis_length::T`: Length of major axis (correlation length)
- `minor_axis_length::T`: Length of minor axis (correlation length)
- `rotation_angle::T`: Rotation angle from geographic east (radians)
- `anisotropy_ratio::T`: Ratio of major to minor axis lengths
"""
struct CorrelationEllipse{T<:AbstractFloat}
    center::Tuple{T,T}
    major_axis_length::T
    minor_axis_length::T
    rotation_angle::T
    anisotropy_ratio::T
    
    function CorrelationEllipse{T}(center::Tuple{T,T}, major_length::T, minor_length::T, 
                                  angle::T) where T
        ratio = major_length / minor_length
        new{T}(center, major_length, minor_length, angle, ratio)
    end
end

"""
    FlowDirectionField{T<:AbstractFloat}

Flow direction analysis results containing directional and magnitude information.

# Fields
- `stream_function::Matrix{T}`: Stream function field ψ
- `flow_angles::Matrix{T}`: Flow direction angles (radians from east)
- `flow_magnitudes::Matrix{T}`: Flow speed magnitudes
- `deformation_tensor::Array{T,3}`: Flow deformation tensor components
- `anisotropy_field::Matrix{T}`: Computed anisotropy ratios
"""
struct FlowDirectionField{T<:AbstractFloat}
    stream_function::Matrix{T}
    flow_angles::Matrix{T}
    flow_magnitudes::Matrix{T}
    deformation_tensor::Array{T,3}
    anisotropy_field::Matrix{T}
end

"""
    StreamFunctionAnalysis{T<:AbstractFloat}

Analyzes atmospheric flow fields to compute stream function and flow characteristics.

# Fields
- `config::AnisotropyConfiguration{T}`: Configuration parameters
- `grid_spacing::Tuple{T,T}`: Grid spacing (dx, dy) in meters
- `laplacian_operator::SparseMatrixCSC{T,Int}`: Laplacian operator for ψ computation
- `gradient_operators::Tuple{SparseMatrixCSC{T,Int}, SparseMatrixCSC{T,Int}}`: ∇x, ∇y operators
"""
struct StreamFunctionAnalysis{T<:AbstractFloat}
    config::AnisotropyConfiguration{T}
    grid_spacing::Tuple{T,T}
    laplacian_operator::SparseMatrixCSC{T,Int}
    gradient_operators::Tuple{SparseMatrixCSC{T,Int}, SparseMatrixCSC{T,Int}}
    
    function StreamFunctionAnalysis{T}(config::AnisotropyConfiguration{T},
                                      grid_size::Tuple{Int,Int},
                                      grid_spacing::Tuple{T,T} = (T(50000), T(50000))) where T
        nx, ny = grid_size
        
        # Create finite difference operators
        lapl_op = create_laplacian_operator(T, nx, ny, grid_spacing)
        grad_x, grad_y = create_gradient_operators(T, nx, ny, grid_spacing)
        
        new{T}(config, grid_spacing, lapl_op, (grad_x, grad_y))
    end
end

"""
    AnisotropicFilter{T<:AbstractFloat}

Main anisotropic filtering operator combining flow analysis and adaptive correlation.

# Fields
- `config::AnisotropyConfiguration{T}`: Configuration parameters
- `stream_analysis::StreamFunctionAnalysis{T}`: Stream function analyzer
- `correlation_ellipses::Matrix{CorrelationEllipse{T}}`: Ellipse field
- `flow_field::Union{FlowDirectionField{T}, Nothing}`: Current flow analysis
- `anisotropy_operators::Dict{Tuple{Int,Int}, SparseMatrixCSC{T,Int}}`: Local operators
"""
mutable struct AnisotropicFilter{T<:AbstractFloat}
    config::AnisotropyConfiguration{T}
    stream_analysis::StreamFunctionAnalysis{T}
    correlation_ellipses::Matrix{CorrelationEllipse{T}}
    flow_field::Union{FlowDirectionField{T}, Nothing}
    anisotropy_operators::Dict{Tuple{Int,Int}, SparseMatrixCSC{T,Int}}
    
    function AnisotropicFilter{T}(config::AnisotropyConfiguration{T},
                                 grid_size::Tuple{Int,Int}) where T
        nx, ny = grid_size
        
        # Initialize stream function analysis
        stream_analyzer = StreamFunctionAnalysis{T}(config, grid_size)
        
        # Initialize with isotropic ellipses
        ellipses = Matrix{CorrelationEllipse{T}}(undef, nx, ny)
        base_length = config.base_correlation_length
        
        for j in 1:ny, i in 1:nx
            center = (T(i), T(j))
            ellipses[i, j] = CorrelationEllipse{T}(center, base_length, base_length, T(0))
        end
        
        operators = Dict{Tuple{Int,Int}, SparseMatrixCSC{T,Int}}()
        
        new{T}(config, stream_analyzer, ellipses, nothing, operators)
    end
end

AnisotropicFilter(config::AnisotropyConfiguration{T}, grid_size::Tuple{Int,Int}) where T = 
    AnisotropicFilter{T}(config, grid_size)

"""
    compute_flow_direction(analyzer::StreamFunctionAnalysis{T}, 
                          u_wind::Matrix{T}, v_wind::Matrix{T}) where T

Compute flow direction field from wind components.

# Arguments
- `analyzer::StreamFunctionAnalysis{T}`: Stream function analyzer
- `u_wind::Matrix{T}`: Zonal wind component (m/s)
- `v_wind::Matrix{T}`: Meridional wind component (m/s)

# Returns
- `FlowDirectionField{T}`: Complete flow analysis including stream function

# Algorithm
1. Compute stream function from vorticity: ∇²ψ = ζ = ∂v/∂x - ∂u/∂y
2. Derive flow angles from stream function gradients
3. Compute flow deformation tensor for anisotropy analysis
4. Determine local anisotropy ratios from deformation characteristics
"""
function compute_flow_direction(analyzer::StreamFunctionAnalysis{T}, 
                               u_wind::Matrix{T}, v_wind::Matrix{T}) where T
    nx, ny = size(u_wind)
    grad_x, grad_y = analyzer.gradient_operators
    
    # Compute vorticity: ζ = ∂v/∂x - ∂u/∂y
    u_flat = vec(u_wind)
    v_flat = vec(v_wind)
    
    dudx = grad_x * u_flat
    dvdy = grad_y * v_flat
    dudy = grad_y * u_flat  
    dvdx = grad_x * v_flat
    
    vorticity = dvdx - dudy
    
    # Solve for stream function: ∇²ψ = ζ
    # Using sparse linear solver (would use iterative methods for large systems)
    psi_flat = analyzer.laplacian_operator \ vorticity
    stream_function = reshape(psi_flat, nx, ny)
    
    # Apply smoothing to stream function
    smooth_factor = analyzer.config.stream_function_filtering
    stream_function = smooth_field(stream_function, smooth_factor)
    
    # Compute flow angles from stream function gradients
    # u = -∂ψ/∂y, v = ∂ψ/∂x
    psi_flat_smooth = vec(stream_function)
    dpsi_dx = grad_x * psi_flat_smooth
    dpsi_dy = grad_y * psi_flat_smooth
    
    u_geostrophic = -reshape(dpsi_dy, nx, ny)
    v_geostrophic = reshape(dpsi_dx, nx, ny)
    
    # Flow angles (from geographic east, counter-clockwise positive)
    flow_angles = atan.(v_geostrophic, u_geostrophic)
    
    # Flow magnitudes
    flow_magnitudes = sqrt.(u_geostrophic.^2 + v_geostrophic.^2)
    
    # Compute deformation tensor for anisotropy analysis
    deformation_tensor = compute_deformation_tensor(grad_x, grad_y, u_wind, v_wind, nx, ny)
    
    # Compute anisotropy field from deformation characteristics
    anisotropy_field = compute_anisotropy_from_deformation(deformation_tensor, 
                                                          analyzer.config.anisotropy_ratio_max)
    
    return FlowDirectionField{T}(stream_function, flow_angles, flow_magnitudes, 
                                deformation_tensor, anisotropy_field)
end

"""
    compute_deformation_tensor(grad_x, grad_y, u_wind, v_wind, nx, ny)

Compute the 2D deformation tensor components for anisotropy analysis.
"""
function compute_deformation_tensor(grad_x, grad_y, u_wind::Matrix{T}, v_wind::Matrix{T}, 
                                   nx::Int, ny::Int) where T
    u_flat = vec(u_wind)
    v_flat = vec(v_wind)
    
    # Velocity gradients
    dudx = reshape(grad_x * u_flat, nx, ny)
    dudy = reshape(grad_y * u_flat, nx, ny)  
    dvdx = reshape(grad_x * v_flat, nx, ny)
    dvdy = reshape(grad_y * v_flat, nx, ny)
    
    # Deformation tensor components
    deformation = zeros(T, nx, ny, 4)
    
    # Stretching deformation: D₁₁ = ∂u/∂x, D₂₂ = ∂v/∂y
    deformation[:, :, 1] = dudx
    deformation[:, :, 2] = dvdy
    
    # Shearing deformation: D₁₂ = D₂₁ = (1/2)(∂u/∂y + ∂v/∂x)
    shear = 0.5 * (dudy + dvdx)
    deformation[:, :, 3] = shear
    deformation[:, :, 4] = shear
    
    return deformation
end

"""
    compute_anisotropy_from_deformation(deformation_tensor::Array{T,3}, max_ratio::T) where T

Compute anisotropy ratios from flow deformation characteristics.
"""
function compute_anisotropy_from_deformation(deformation_tensor::Array{T,3}, max_ratio::T) where T
    nx, ny, _ = size(deformation_tensor)
    anisotropy = ones(T, nx, ny)
    
    for j in 1:ny, i in 1:nx
        # Extract 2x2 deformation tensor
        D = @SMatrix [deformation_tensor[i,j,1] deformation_tensor[i,j,3];
                      deformation_tensor[i,j,4] deformation_tensor[i,j,2]]
        
        # Compute eigenvalues to determine principal deformation directions
        eigenvals = eigvals(D)
        λ₁, λ₂ = real(eigenvals[1]), real(eigenvals[2])
        
        # Anisotropy ratio based on eigenvalue ratio
        if abs(λ₂) > 1e-10
            ratio = abs(λ₁ / λ₂)
            anisotropy[i, j] = min(ratio, max_ratio)
        else
            anisotropy[i, j] = T(1)  # Isotropic if no deformation
        end
    end
    
    return anisotropy
end

"""
    update_correlation_ellipses!(filter::AnisotropicFilter{T}, 
                                flow_field::FlowDirectionField{T}) where T

Update correlation ellipses based on flow field analysis.
"""
function update_correlation_ellipses!(filter::AnisotropicFilter{T}, 
                                     flow_field::FlowDirectionField{T}) where T
    nx, ny = size(flow_field.flow_angles)
    base_length = filter.config.base_correlation_length
    adaptation_factor = filter.config.flow_adaptation_factor
    
    for j in 1:ny, i in 1:nx
        # Flow-adapted correlation lengths
        flow_magnitude = flow_field.flow_magnitudes[i, j]
        anisotropy_ratio = flow_field.anisotropy_field[i, j]
        
        # Major axis aligned with flow direction
        flow_angle = flow_field.flow_angles[i, j]
        
        # Adapt correlation lengths based on flow characteristics
        major_length = base_length * (1 + adaptation_factor * anisotropy_ratio)
        minor_length = base_length / sqrt(anisotropy_ratio)
        
        # Apply rotation smoothing
        if i > 1 && j > 1
            prev_angle = filter.correlation_ellipses[i-1, j].rotation_angle
            smooth_factor = filter.config.ellipse_rotation_smoothing
            smoothed_angle = (1 - smooth_factor) * flow_angle + smooth_factor * prev_angle
        else
            smoothed_angle = flow_angle
        end
        
        # Update ellipse
        center = (T(i), T(j))
        filter.correlation_ellipses[i, j] = CorrelationEllipse{T}(center, major_length, 
                                                                 minor_length, smoothed_angle)
    end
    
    # Store current flow field
    filter.flow_field = flow_field
end

"""
    apply_anisotropic_filter(filter::AnisotropicFilter{T}, field::Matrix{T}, 
                            u_wind::Matrix{T}, v_wind::Matrix{T}) where T

Apply anisotropic filtering with flow-dependent correlation structure.

# Arguments
- `filter::AnisotropicFilter{T}`: Configured anisotropic filter
- `field::Matrix{T}`: Input field to be filtered  
- `u_wind::Matrix{T}`: Zonal wind component for flow analysis
- `v_wind::Matrix{T}`: Meridional wind component for flow analysis

# Returns
- `Matrix{T}`: Filtered field with anisotropic correlations

# Algorithm
1. Analyze flow field to determine correlation ellipse parameters
2. Update correlation ellipses based on flow characteristics
3. Apply directional filtering aligned with correlation ellipses
4. Combine multi-directional filter results with proper weighting
"""
function apply_anisotropic_filter(filter::AnisotropicFilter{T}, field::Matrix{T}, 
                                 u_wind::Matrix{T}, v_wind::Matrix{T}) where T
    # Analyze flow field
    flow_field = compute_flow_direction(filter.stream_analysis, u_wind, v_wind)
    
    # Update correlation ellipses
    update_correlation_ellipses!(filter, flow_field)
    
    # Apply anisotropic filtering
    result = apply_ellipse_based_filtering(filter, field)
    
    return result
end

"""
    apply_ellipse_based_filtering(filter::AnisotropicFilter{T}, field::Matrix{T}) where T

Apply filtering based on correlation ellipse field.
"""
function apply_ellipse_based_filtering(filter::AnisotropicFilter{T}, field::Matrix{T}) where T
    nx, ny = size(field)
    result = copy(field)
    temp_result = similar(field)
    
    # Apply anisotropic filtering at each grid point
    for j in 1:ny, i in 1:nx
        ellipse = filter.correlation_ellipses[i, j]
        
        # Create local anisotropic operator if not cached
        ellipse_key = (i, j)
        if !haskey(filter.anisotropy_operators, ellipse_key)
            filter.anisotropy_operators[ellipse_key] = 
                create_local_anisotropic_operator(ellipse, nx, ny, i, j)
        end
        
        # Apply local anisotropic filtering
        local_operator = filter.anisotropy_operators[ellipse_key]
        field_flat = vec(field)
        filtered_flat = local_operator * field_flat
        temp_result .= reshape(filtered_flat, nx, ny)
        
        # Weight contribution based on distance from grid point
        weight = compute_local_weight(ellipse, i, j, nx, ny)
        result[i, j] = weight * temp_result[i, j] + (1 - weight) * field[i, j]
    end
    
    return result
end

"""
    create_local_anisotropic_operator(ellipse::CorrelationEllipse{T}, nx::Int, ny::Int, 
                                     center_i::Int, center_j::Int) where T

Create sparse operator for local anisotropic filtering around an ellipse.
"""
function create_local_anisotropic_operator(ellipse::CorrelationEllipse{T}, nx::Int, ny::Int,
                                          center_i::Int, center_j::Int) where T
    n_total = nx * ny
    
    # Create sparse matrix for local anisotropic operator
    I_indices = Int[]
    J_indices = Int[]  
    values = T[]
    
    # Rotation matrix for ellipse orientation
    θ = ellipse.rotation_angle
    cos_θ, sin_θ = cos(θ), sin(θ)
    rotation_matrix = @SMatrix [cos_θ -sin_θ; sin_θ cos_θ]
    
    # Correlation length scales
    a = ellipse.major_axis_length  # Major axis length
    b = ellipse.minor_axis_length  # Minor axis length
    
    # Define influence radius (3 times major axis length)
    influence_radius = 3 * a
    
    for j in 1:ny, i in 1:nx
        target_idx = (j-1)*nx + i
        
        # Compute distance from ellipse center
        dx = T(i - center_i)
        dy = T(j - center_j)
        
        # Skip if outside influence radius
        if sqrt(dx^2 + dy^2) > influence_radius
            continue
        end
        
        # Rotate to ellipse coordinate system
        rotated_coords = rotation_matrix * @SVector [dx, dy]
        x_rot, y_rot = rotated_coords[1], rotated_coords[2]
        
        # Compute anisotropic correlation weight
        # Using elliptical Gaussian: exp(-(x²/a² + y²/b²)/2)
        ellipse_distance_squared = (x_rot/a)^2 + (y_rot/b)^2
        weight = exp(-ellipse_distance_squared / 2)
        
        # Add to sparse matrix if weight is significant
        if weight > 1e-6
            push!(I_indices, target_idx)
            push!(J_indices, target_idx)
            push!(values, weight)
        end
    end
    
    # Normalize weights to preserve variance
    total_weight = sum(values)
    if total_weight > 0
        values ./= total_weight
    end
    
    return sparse(I_indices, J_indices, values, n_total, n_total)
end

"""
    compute_local_weight(ellipse::CorrelationEllipse{T}, i::Int, j::Int, 
                        nx::Int, ny::Int) where T

Compute local weighting for ellipse-based filtering.
"""
function compute_local_weight(ellipse::CorrelationEllipse{T}, i::Int, j::Int, 
                             nx::Int, ny::Int) where T
    # Simple distance-based weighting - could be made more sophisticated
    center_x, center_y = ellipse.center
    distance = sqrt((T(i) - center_x)^2 + (T(j) - center_y)^2)
    max_influence = ellipse.major_axis_length
    
    if distance > max_influence
        return T(0)
    else
        return exp(-distance^2 / (2 * max_influence^2))
    end
end

"""
    smooth_field(field::Matrix{T}, smoothing_factor::T) where T

Apply simple smoothing to a field.
"""
function smooth_field(field::Matrix{T}, smoothing_factor::T) where T
    if smoothing_factor <= 0
        return field
    end
    
    nx, ny = size(field)
    smoothed = copy(field)
    
    # Simple 3x3 kernel smoothing
    for j in 2:(ny-1), i in 2:(nx-1)
        neighborhood_mean = (field[i-1,j-1] + field[i,j-1] + field[i+1,j-1] +
                           field[i-1,j] + field[i,j] + field[i+1,j] +
                           field[i-1,j+1] + field[i,j+1] + field[i+1,j+1]) / 9
        
        smoothed[i,j] = (1 - smoothing_factor) * field[i,j] + 
                        smoothing_factor * neighborhood_mean
    end
    
    return smoothed
end

"""
    create_laplacian_operator(::Type{T}, nx::Int, ny::Int, 
                             grid_spacing::Tuple{T,T}) where T

Create sparse Laplacian operator for 2D grid.
"""
function create_laplacian_operator(::Type{T}, nx::Int, ny::Int, 
                                  grid_spacing::Tuple{T,T}) where T
    dx, dy = grid_spacing
    n = nx * ny
    
    I_indices = Int[]
    J_indices = Int[]
    values = T[]
    
    for j in 1:ny, i in 1:nx
        idx = (j-1)*nx + i
        
        # Central differences for interior points
        if 1 < i < nx && 1 < j < ny
            # Central point
            push!(I_indices, idx); push!(J_indices, idx); push!(values, -2/dx^2 - 2/dy^2)
            
            # x-neighbors
            push!(I_indices, idx); push!(J_indices, idx-1); push!(values, 1/dx^2)
            push!(I_indices, idx); push!(J_indices, idx+1); push!(values, 1/dx^2)
            
            # y-neighbors  
            push!(I_indices, idx); push!(J_indices, idx-nx); push!(values, 1/dy^2)
            push!(I_indices, idx); push!(J_indices, idx+nx); push!(values, 1/dy^2)
        else
            # Boundary points - simplified zero gradient
            push!(I_indices, idx); push!(J_indices, idx); push!(values, T(1))
        end
    end
    
    return sparse(I_indices, J_indices, values, n, n)
end

"""
    create_gradient_operators(::Type{T}, nx::Int, ny::Int, 
                             grid_spacing::Tuple{T,T}) where T

Create sparse gradient operators for 2D grid.
"""
function create_gradient_operators(::Type{T}, nx::Int, ny::Int, 
                                  grid_spacing::Tuple{T,T}) where T
    dx, dy = grid_spacing
    n = nx * ny
    
    # X-gradient operator
    I_x, J_x, V_x = Int[], Int[], T[]
    # Y-gradient operator  
    I_y, J_y, V_y = Int[], Int[], T[]
    
    for j in 1:ny, i in 1:nx
        idx = (j-1)*nx + i
        
        # X-gradient
        if i > 1 && i < nx
            # Central difference
            push!(I_x, idx); push!(J_x, idx-1); push!(V_x, -T(0.5)/dx)
            push!(I_x, idx); push!(J_x, idx+1); push!(V_x, T(0.5)/dx)
        elseif i == 1
            # Forward difference
            push!(I_x, idx); push!(J_x, idx); push!(V_x, -T(1)/dx)
            push!(I_x, idx); push!(J_x, idx+1); push!(V_x, T(1)/dx)
        else # i == nx
            # Backward difference
            push!(I_x, idx); push!(J_x, idx-1); push!(V_x, -T(1)/dx)
            push!(I_x, idx); push!(J_x, idx); push!(V_x, T(1)/dx)
        end
        
        # Y-gradient
        if j > 1 && j < ny
            # Central difference
            push!(I_y, idx); push!(J_y, idx-nx); push!(V_y, -T(0.5)/dy)
            push!(I_y, idx); push!(J_y, idx+nx); push!(V_y, T(0.5)/dy)
        elseif j == 1
            # Forward difference
            push!(I_y, idx); push!(J_y, idx); push!(V_y, -T(1)/dy)
            push!(I_y, idx); push!(J_y, idx+nx); push!(V_y, T(1)/dy)
        else # j == ny
            # Backward difference
            push!(I_y, idx); push!(J_y, idx-nx); push!(V_y, -T(1)/dy)
            push!(I_y, idx); push!(J_y, idx); push!(V_y, T(1)/dy)
        end
    end
    
    grad_x = sparse(I_x, J_x, V_x, n, n)
    grad_y = sparse(I_y, J_y, V_y, n, n)
    
    return grad_x, grad_y
end

"""
    validate_anisotropic_response(filter::AnisotropicFilter{T}, 
                                 test_wind::Tuple{Matrix{T}, Matrix{T}}) where T

Validate anisotropic filter response with test wind field.
Returns analysis of ellipse alignment and correlation structure.
"""
function validate_anisotropic_response(filter::AnisotropicFilter{T}, 
                                      test_wind::Tuple{Matrix{T}, Matrix{T}}) where T
    u_wind, v_wind = test_wind
    nx, ny = size(u_wind)
    
    # Create test impulse
    impulse = zeros(T, nx, ny)
    center_x, center_y = nx÷2, ny÷2
    impulse[center_x, center_y] = T(1)
    
    # Apply anisotropic filter
    response = apply_anisotropic_filter(filter, impulse, u_wind, v_wind)
    
    # Analyze response characteristics
    max_response = maximum(response)
    response_center = findmax(response)[2]
    
    # Compute correlation lengths in different directions
    x_profile = response[center_x, :]
    y_profile = response[:, center_y]
    
    # Find half-maximum points for correlation length estimation
    half_max = max_response / 2
    
    return (response=response, max_value=max_response, center=response_center,
            x_profile=x_profile, y_profile=y_profile, half_max=half_max)
end

end # module AnisotropicFiltering